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ABSTRACT 

Since the violent relaxation in hierarchical merging is incomplete, elliptical galaxies retain a wealth 
of information about their formation pathways in their present-day orbital structure. Recent advances 
in integral field spectroscopy, multi-slit infrared spectroscopy, and triaxial dynamical modeling tech- 
niques have greatly improved our ability to harvest this information. A variety of observational and 
theoretical evidence indicates that gas-rich major mergers play an important role in the formation of 
elliptical galaxies. We simulate 1:1 disk mergers at seven different initial gas fractions (f ga s) ranging 
from to 40%, using a version of the TreeSPH code Gadget-2 that includes radiative heating and cool- 
ing, star formation, and feedback from supernovae and active galactic nuclei. We classify the stellar 
orbits in each remnant and construct radial profiles of the orbital content, intrinsic shape, and ori- 
entation. The dissipationless remnants are typically prolate-triaxial, dominated by box orbits within 
r c ~ 1.5i? e , and by tube orbits in their outer parts. As f gas increases, the box orbits within r c are 
increasingly replaced by a population of short axis tubes (z— tubes) with near zero net rotation, and 
the remnants become progressively more oblate and round. The long axis tube (x— tube) orbits are 
highly streaming and relatively insensitive to f gas , implying that their angular momentum is retained 
from the dynamically cold initial conditions. Outside r c , the orbital structure is essentially unchanged 
by the gas. For f gas > 15%, gas that retains its angular momentum during the merger re-forms a 
disk, that appears in the remnants as a highly streaming 2— tube population superimposed on the hot 
z— tube distribution formed by the old stars. In the 15-20% gas remnants, this population appears as a 
kinematically distinct core (KDC) within a system that is slowly rotating or dominated by minor-axis 
rotation. These remnants show an interesting r esemblance, in both their y elocity maps and intrinsic 
orbital structure, to the KDC galaxy NGC4365 (|van den Bosch et al.|[2008h . At 30-40% gas, the rem- 
nants are rapidly rotating, with sharp embedded disks on ~ lR e scales. We predict a characteristic, 
physically intuitive orbital structure for 1:1 disk merger remnants, with a distinct transition between 1 
and 3i? e that will be readily observable with combined data from the 2D kinematics surveys SAURON 
and SMEAGOL. Our results illustrate the power of direct comparisons between N— body simulations 
and dynamical models of observed systems to constrain theories of galaxy formation. 
Subject headings: stellar dynamics - methods: n-body simulations - galaxies: elliptical and lenticu- 
lar, cD - galaxies: formation - galaxies: interactions - galaxies: kinematics and 
dynamics 



1. INTRODUCTION 

1.1. Motivation 

In the standard ACDM c oncordance cosmology 
(lOstriker fe Steinhardtl 119951: iDodelson et~aH 119961: 
ISpergel et al.1 120071 ) . structure in the Universe grows 
hierarchically, through a progression of smaller bodies 
accreti ng material and merg ing to form larger systems 
(e.g. I White fc Reesl H971 . Cosmological N-body 
simulations starting from a Gaussian random field 
of linear density flu ctuations at high redshift (e.g. 
iSpringel et al.l l2005b| ). togeth er with analytic m odels 
of gravitational collapse (e.g. iBertschingerl [1985?) . have 
given us a fairly clear picture of how dark matter (DM) 
assemble s in the Universe. Exten d ed Press-Sc h echter 
theory (|Press fc Schechterl 119741: iBond et al.l 119911: 
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lLacev fc Cold [L993h provides an analytic formalism for 
computing halo merger rates and assembly histories, 
which matches TV-body simulations r emarkably well 
(|Lacev fc ColdH99l : iGenel et al.ll2009bD . and has been 
used in a var iety of semianalytic models of structure for- 
matio n (e.g. iCole et al.ll2000b iManriaue fc Salvador-Sole 
1996 ; iKauffmann et al.l 119991 : iSomerville fc Primack 



19991 iCroton et al.l 120061: ISomerville et al.l 120081: 



Ste wart et al. 2009b). Simulated DM halos also appear 



to share a universal internal morphology, with density 
and velocity anisotropy profiles similar to the generic 
outcome of violent relaxation fol lowing dissipationless 
collap s e or strong tid al shoc ki ng (IDubinski fc Carlbergj 
199lL _ Havarro et al.l 



19971: IBullock et all I20ln1lbf: 



Navarro et al.l l2008t Miller fc Smith! 119791: Ivan Albadal 



1982t iMcGlvnnl 119841 119901: ISpergel 1 fc Hernauistl 11992 
Huss et al.lll999HMacMillan et al.ll2006UBellovarv et al.l 
2008). 

There is no such simple theory of hierarchical galaxy 
formation, because the luminous components of galax- 
ies are formed through complex baryonic physics. For 
instance radiative heating and cooling, star formation 
(SF) and gas expulsion through stellar winds, energy 
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and momentum feedback from supernovae and active 
galactic nuclei (AGN), ram pressure stripping, and res- 
onant effects in dynamically cold systems, are all im- 
portant ingredients in de termining the baryo n ic struc- 



ture of galaxies (e . g. iCroton et al.l 120061 : iSpringe 



T] EL 

20001 iSpringel et all l2005a| iBest et al.l 120071: ICox et al 
2006U iCeverino fe Klvpml 120091: iMartig fe Bournauc 



2009TlD , Onghia et al.l 120091 ) . Cosmological models of 



galaxy formation must therefore be calibrated with high- 
resolution simulations on galactic scales that study these 
effects in isolation, coupled with direct constraints from 
observations. 

Violent relaxation in mergers tends to drive galaxies 
toward a universal , fully mixed structure (lAlladinl 119651: 
Lvnden-Belll fl967t ISver fe Whitel fl998t iThomas et atl 



20091 ). while gas accretion and stell ar outflows produce 



new dynamically cold components d Rix & Whitel [1990 1: 
Block et al.1 120021: iGenel et all l2009al iDekel et aLlbooa: 



Bournaud fc Elmc greenll2009l; IMartig fc Bournaudll2009l: 
Daddi et al.l 12009). From the hot stellar distributions 



of elliptical galaxies, we deduce that they are the most 
merger-dominated systems. Ellipticals and bulges con- 
tain the m ajority of the stellar mass in the local Uni- 
verse (e.g. lGadottil[2009l ). and often serve as a testing 
ground for theory since they are the most evolved under 
the complex combination of processes driving galaxy for- 
mation. Since the violent relaxation in mergers is incom- 
plete, elliptical galaxies retain a wealth of information 
about their formation hi stories in their present-day dis- 
tribution functions (e.g. lEggen et aT1ll962l : iLvnden-Belll 
IT967l:IWhitelll980HValluri et alJl2007D . 

Recent advances in integral field spectroscopy (IFS : 
Bacon et al.lU995. 2001:lHill et"aI|[2006l:IWeiimans et al 



2009), multi-slit inf r ared spectroscopy ()Norris et al 
20081: IProctor et all 120091 ), and triaxial d y nami- 
cal modeling techniques (Ivan de Ven et al.l 1200 



Ivan den Bosch et al.l 120081: Ide Lorenzi et al.l 120061 I200S 
have greatly improved our ability to harvest this in- 
formation. The SAU RON project (|Bacon et al.l 120011: 
lEmsellem et al.l 120041 ) will produce high resolution, 
2D kinematic maps within ~ 1 effective radius (R e ) 
for a representative sample of ^100 nearby elliptical 
galaxies and spiral bulges, using a panoramic integral 
field spectrograph mounted on the William Herschel 
Telescope. The data on 48 early-type galaxies released 
to date has revealed an unexpectedly rich variety of 
kinematic structures, which poses a new challenge 
for galaxy formati on simulations (Jesseit e t al.l 120071 : 
iBurkert et al.l [20081. Dynamical modeling studies have 
shown that, in practice, 2D maps of the first four mo- 
ments (hi — h^) of the line-of-sight velocity distribution 
(LOSVD) provided by SAURON are typically sufficient 
to uniqu ely reconstruct t h e 3D s tellar orbital distri- 
bution (van de Ven et all l200l Ivan den Bosch et al.l 
|2008[ ). Complex features present in many systems, such 
as embedded disks and kinematically distinct cores 
(KDCs), can provi de especially strong constraints on th e 
intrinsic structure (jvan den Bosch fc van de Venll2008l ). 

A good example of the power of these new techniques is 
provided by the case of NGC4365. This m assive old ellip- 
tical is known for its min or axis rot ation dWagner et al.l 
fl988tlBender et al.lll99l and KDC (IDavies et al.ll2001l) . 
and is therefore a na t ural c andidate for dynamical mod- 
eling. iStatler et al.l (|2004D modeled this galaxy using 



a velocity field fitting method (|Statlerl 1994a, b) that 
made use of the surface brightness and full 2D velocity 
map from SAURON, but not the higher moments of the 
LOSVDs. They found that the system was nearly max- 
imally triaxial, ruling out axisymmetry at >95% confi- 
dence. 

A few years later, Ivan den Bosch etlrtl (|2008l) mod- 
eled the same galaxy u sing an advanced n ew triaxial 
Schwarzschild modeling (Schwarzschild 1979) code that 
can incorporate all of the LOSVD moments up to /14. 
They reached a qualitatively different conclusion - that 
the system was nearly oblate axisymmetric. The pre- 
dominance of the minor-axis rotation in the outer parts 
of the map owed to a high degree of cancellation of the 
orbits rotating in a prograde and retrograde sense about 
the short axis, and streaming of the smaller population 
of orbits rotating about the intrinsic long axis. They 
found no major transition in the orbital structure at the 
boundary of the kinematically "decoupled" core, mak- 
ing it unlikely that it formed in a separate infall event. 
The orbital structure of the 15-20% gas merger remnants 
in this paper bears a tantalizing resemblance to that of 
NGC4365 (compare e.g. ou r Figures ®1 - ITDI and IT71 - ITgl 
with Figures 7, 11 and 12 of Ivan den Bosch et al. 2008). 

However a limitation of the SAURON spectrograph 
is its small field of view, corresponding to ~ lR e on a 
typical elliptical target. The outer parts of galaxies are 
less relaxed than their inner parts, and bar-like modes 
in mer gers efficiently transport angula r mom e ntum out 



1973 



1996; 



war d (IQstriker fc Peebles! 
1993: iBarnes fc Hcrn quistl 

Hoff man et al.l l2009cf) . A gas- rich merger between two 



Fall [19791: IHernquist] 
Hop kins et alT l2008b: 



spiral galaxies with halos might be expected to produce 
a remnant with three distinct components in its distri- 
bution function: (i) an inner part formed through dissi- 
pation, (ii) a middle part reflecting the dynamically cold 
distribution of the disk stars, and (iii) an outer part aris- 
ing from the pre-existing stellar halo populations. Ob- 
servations with about four times the spatial coverage of 
SAURON would be able to detect these dynamical sub- 
components, and probe the parts of galaxies retaining the 
most memory of their progenitors' angular momentum 
and internal structure. Hints of increased complexity in 
the angular momentum profiles at large ra dii have indeed 
been observed in a few el liptical galaxies (jCoccato et al.l 
f2009t IProctor et al.ll200l . 

A number of projects designed to extend SAURON- 
style dynamical modeling out to larger radii are cur- 
rently underway. The S MEAGOL survey ()Proctor et al.l 
I2009t iFoster et al.ll2009l) will obtain smoothed 2D maps 
of hi — /14 out to ~ 3i? e for a representative sample 
of 25 nearby ellipticals, using the new stel lar kinemat- 
ics w i th multiple slits (SK i MS) t echnique (jNorris et al.l 
120081 : IProctor et ail 120081 f2009h with the DEIMOS 
spectrograph on the 10-meter Keck-II telescope. The 
data analysis will include triaxial dynamical model- 
ing with the advanced p a rticle-based method NMAGIC 
(ISver fc Tremaind fl996i: Ide Lorenzi etaLl [200l I2008L 
The wide-field IFS VIRUS-P has been used 



to obtain 2D stellar kin ematics out to 3 — AR P for 
severa l giant ellipticals (jHill et all 120061 : iBlanc et~al"1 
I2009t iMurphv et all |2009|), and multiple pointings of 
the SAURON spectrograph have been used to measure 
the stellar LOSVDs out to ~ 4i? e in NGC3379 and 
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NGC821 (IWeiimans et al.ll2009fl . Surveys using globular 
clusters (GCs) and planetary nebul ae (PNe) as discrete 
tracers of the distribution funct i on (iRomanowskv et al.l 
2009aHDouglas et al.1120021 [20071: iNapolitano et al ll2009t 
Nantais fc Huchral 120091 ISchuberth et al.l I2009Q can 
probe much larger radii (out to ~ 10 R e ) and place fur- 
ther constraints on dynamical models (Chanam e et al.l 
l2008h . 

Combined with simulations aimed at establishing the 
characteristic orbital structure arising from various for- 
mation pathways, these observational programs will pro- 
vide unprecedented insight into the physics of galaxy for- 
mation. 

1.2. Orbits in galactic potentials 

To deduce the formation histories of galaxies from their 
orbital structure, we must identify macroscopic groups 
of orbits that are confined to a well-defined neighbor- 
hood of phase space because they followed similar evo- 
lutionary pathways. The isolating integrals of motion 
(or quas i-isolating integrals i n the case of perturbed po- 
tentials; IContopoul os 1963b; jGoodman fc Schwarzschildl 
[l98lt ISiopis fc Kandrupl I2000t IKandrup fc Siopisl I200"3¥ 
parameterize the phase space region in which a star re- 
mains localized, and encode whatever information about 
its initial conditions (ICs) is preserved once the sys- 
tem is fully phase-mixed (e.g.. IBinnev fc Tremainel 120081 : 
IBinnev fc Sp"ergel [l98l IGomez fc Helmil 120091) . Orbits 
which conserve at least one isolating integral per degree 
of freedom are called regular. It is the ubiquity of reg- 
ular orbits in galaxies that permits the rich variety in 
their structure, mirroring their varied formation path- 
ways (jSchwarzschildlll979t) . 

The regular orbits in a static potential can be classified 
into families that conserve qualitatively similar integrals 
of the motion, and therefore have similar morphologies. 
Which orbital class a given star will occupy is determined 
by the available phase space for different kinds of orbits 
in the potential, and its ICs - the phase space need not 
be uniformly populated. In a time-varying (or otherwise 
non-ideal) potential such as an ongoing merger, stars dif- 
fuse in the space of their conserved integrals, but not 
completely. If they do not cross boundaries between the 
orbital families, then their qualitative character may be 
preserved from the ICs. Crossing between orbital bound- 
aries can serve as a collect ive relaxation mechanism (e.g. 
IBarnes fc Hern quist 1991) ■ An intuitive grasp of the or- 
bital classes is therefore essential to understanding how 
dynamical systems evolve and relax. 

We begin with a brief overview of the types of regular 
orbits that are possible in various idealized potentials, 
leading up to a classification of the orbits in triaxial sys- 
tems into families that conserve similar integrals. For 
a more t horough and rigorous prese ntation we refer the 
reader to IBinnev fc Tremainel (|2008l ). 

The simplest conceivable model is a spherical poten- 
tial. In this case the symmetry about all three Carte- 
sian axes implies conservation of the angular momentum 
vector, so every orbit is confined to a plane. The star 
oscillates in radius with frequency w r while precessing in 
azimuth with frequency ujg. If these two frequencies are 
commensurate (mw r -I- nuig = for some integers m and 
n) then the orbit closes on itself, as in a Kepler potential. 
More generally, the orbits form rosettes that eventually 



fill an annulus between the minimum and maximum of 
the radial oscillations (pericenter and apocenter). 

Very few galaxies are spherical, but many are consis- 
tent with axisymmetry. In an axisymmetric potential 
the angular momentum component about the symme- 
try axis, L Zy is conserved. The direction of L precesses 
about the z axis as L x and L y vary. The radial oscilla- 
tions no longer return the star to the exact same peri- 
center and apocenter every cycle, but are still bounded 
between some p = r m i n and a — r max . 

Two orbits with the same energy and L z can look quite 
different from each other, ranging from orbits nearly con- 
fined to the x — y plane, resembling eccentric orbits in 
a thin disk, to puffed-up orbits nearly filling a spherical 
annulus over long times. This suggests that the orbits 
are constrained by another integral in addition to E and 
L z , related to how the energy is apportioned into vertical 
and radial motion. Though it cannot be expressed ana- 
lytically for a general axisymmetric potential, this third 
integral (1$) allows an assortment of stable axisymmetric 
systems, from thin disks formed by quiescent accretion, 
to nearly spherical systems heated b y persistent pertur- 
bations or discrete encounters (e.E 



1963aHHenon fc Heileslll964l; lSaaflll968l; lRichstonelll98 



Bin nev fc Tremainel l2008) . 
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Fig. 1. — Toy model of the box-tube boundary: Elliptical 
ring in a frictionless, concave trough. The ring is initially 
in equilibrium, with its long axis aligned with the long axis of the 
trough. If an angular velocity lj is imparted to the ring, it must 
go up and over a barrier to spin all the way around. Above some 
critical uj c , the ring will make it over the barrier and continue spin- 
ning; below oj c it will just librate about the eq uilibrium position. 
This analogue owes to Binnev & Spcrgcl (1982). 

However some galaxies show clear evidenc e of nonax- 
isymmetric shapes (e.g. iFranx et al.l Il99l[) . and sim- 
ulations of dissipationless violent relaxation in merg- 
ers or collapses generica lly produce triaxial remnants 
(e.g. Ivan Albadal fl98l . There is a unique den- 
sity distribution stratified on similar triaxial ellipsoids 
whose potential is separable in ellipsoidal coordinates 
(|de Zeeuw fc Lvnden-Belll |l985), known as the "per- 
fect ellipsoid". Its distribution function and integrals 
of motion can be expressed analytically, and its or- 
bital structure has be e n studied extensively dStackell 
T8901: lEddinetonl H9l5t fde Zeeuw fc Lvnden-Belll !T98a 



de Zeeuwl 119851 : iStatlerl fl987t ). The orbits in the per 
feet ellips o id we re classified into four major families by 
Ide Zeeuwl (|1985f ): short-axis tubes (z— tubes) which ro- 
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tate about the short (z— )axis of the potential, two classes 
of long-axis tubes (x— tubes) which revolve about the 
long (x— )axis, and box orbits which behave like per- 
turbed simple harmonic oscillators. Because it is ana- 
lytic, this model has been widely used as a starting point 
for understanding the orbital structure of more general 
triaxial potentials. 

Tube orbits resemble the orbits in axisymmetric po- 
tentials, and may be thought of loosely as precessing 
ellipses drive n by the bar-like poten tial of the triax- 
ial ellipsoid (jBinnev & Spergell 11983) . They conserve 
angular momentum-like integrals and therefore avoid 
the origin of the potential and the zero-velocity sur- 
face. It can be shown that y— tube orbits rotating about 
the intermediate axis in the perfect ellipsoid are unsta- 
ble to vertical perturba tions ( Hciligman fc Schwarzschildl 
Il979t Ide Zeeuwj[i985D . but both x-tube and z-tube 
orbits are allowed. Observations of both major- and 
minor-axis rotation in some elliptical galaxies there- 
fore strongly suggests that these systems are triaxial 
(IIllingwortrj'(l977l:lSchechter fc Gunnlll979t lBilmevlll985l: 
IWagner et all 119881: iFranx et al.l I1991H T x— tube orbits 
are most prevalent in prolate potentials, and are popu- 
lated by stars with large initial angular momenta about 
the long axis, z— tubes are the dominant type of orbit 
in oblate systems. Tube orbits oscillate in radius within 
some bounds p and a, and one component of L never 
switches sign. Any net angular momentum of a triax- 
ial system must be carried by the tube orbits, so only 
these orbits can retain information about a galaxy's ini- 
tial sense of rotation. 

Just as tube orbits may be thought of as precessing 
ellipses, box orbits may be regarded loosely as axial or- 
bits (or elongated ellips es) librating about the a;— axis. 
iBinnev fc Spergel (|1982f ) illustrated the transition be- 
tween box and tube orbits using an intuitive toy model. 
Imagine an elliptical ring lying in a concave, frictionless 
trough, with its long axis initially parallel to the long axis 
of the trough, as shown in Figure [TJ This configuration 
allows the ring to lie as low in the trough as possible. 
Now imagine trying to spin the ring in the trough, by 
applying an impulsive kick of energy \lw 2 to one of its 
ends. To spin all the way around the ring must go up and 
over a barrier, since the curvature of the trough along the 
perpendicular direction is greater. There is some criti- 
cal precession frequency u) c , above which the ring will 
get over the barrier. Below uj c it will just librate about 
the equilibrium configuration, with no definite sense of 
rotation. 

The trough's different curvature along the two axes is 
analogous to triaxiality of a gravitational potential, and 
the librating mode is analogous to box orbits. If the cur- 
vature is the same along both directions perpendicular 
to the rotation axis (as in an axisymmetric potential), 
then the energy barrier is zero. Note also that the effec- 
tive barrier in a triaxial potential is greatest for a star 
rotating about the y— axis, since in this case the differ- 
ence in curvature along the two perpendicular directions 
is greatest, making y— tubes the most susceptible to in- 
stability. 

Box orbits are prevalent in triaxial systems with shal- 
low inner density profiles, and conserve integrals simi- 
lar to the energies of independent harmonic oscillations 



about each Cartesian axis. Stars on box orbits have no 
definite sense of rotation, and can therefore pass arbi- 
trarily close to both the origin (they are "centrophilic" ) 
and the zero-velocity surface. Over time, they densely 
fill a 3-dimensional box-like region centered on the origi n 
(|de Zeeuwlll985t IStatlerill987t IBinnev fc Tremainef2 008) . 
Powerful nuclear processes such as gas inflow, starbursts, 
and black hole growth are thought to be major drivers of 
galaxy evolution, and box orbits may be responsible for 
conveying information about the rapidly varying central 
potential to large radii. 
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Fig. 2. — Orbital composition of a typical merger rem- 
nant. Remnant i (see Table 1), 15% gas, viewed along a sight-line 
near the y— axis. First column: Projected surface brightness, E, 
colored on a logarithmic scale. All four maps use the same color 
scaling to show the relative contribution of each orbital class. Sec- 
ond column: LOS velocity v, with overplotted isophotes. The size 
of the maps is lR e X lR e , and the color scale runs from —v max to 
Umax, where v max is the maximum value of |i;| in any pixel. The 
number in the upper right corner gives the value of Vmax in km/s. 
The two right panels show the smoothed LOSVDs at the locations 
specified by the white letters on the E and v maps. Top: LOSVD 
at point A, near the major axis. Bottom: LOSVD at point B, near 
the minor axis. The heavy black line is the full LOSVD, while the 
thin lines with markers break it down by orbital class. 

Figure [2] shows an example of how the three orbital 
classes contribute to the projected surface brightness and 
kinematics in one of our simulated merger remnants. The 
remnant is shown in projection along the y— axis, which 
maximally separates the orbital classes in 2D space. The 
boxes and z— tubes are elongated along the major axis, 
while the x— tubes are elongated along the minor axis 
of the projection. All four surface brightness maps are 
plotted on the same color scale, to show the relative or- 
bital populations at different locations in the sky plane. 
Note that the isophotes of each individual orbital class 
appear boxier than the combined isophotes, which are 
nearly elliptical in shape. 

The a;— tube orbits are responsible for the minor-axis 
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rotation in the remnant, and the z— tubes produce the 
major-axis rotation. There is a striking difference in the 
amount of streaming of the two classes of tubes: the 
velocity scale on the a:— tube map is ±145 km/s, while 
it is only ±53 km/s on the z— tube map. The remnant 
therefore appears dominated by minor-axis rotation in 
the velocity maps, even though z— tube orbits dominate 
its stellar mass. 

Two local LOSVDs are shown on the right, at the loca- 
tions indicated on the maps: one at a point near the cen- 
ter along the major axis (point A) and one near the minor 
axis farther out (point B). At point A the z— tubes have 
a flat-topped distribution with a high velocity dispersion 
and negative kurtosis, indicative of canceling streams ro- 
tating in opposite directions. The box orbits are strongly 
peaked at v — 0, with a large excess in the tails giving 
the distribution a high positive kurtosis. The combi- 
nation of the z— tube and box orbits yields a combined 
LOSVD that is closer to a Gaussian. At point B the 
a;— tube orbits make a larger relative contribution. Their 
LOSVD is peaked at a high negative velocity and skewed 
right, indicative of a heated streaming population. The 
combined distribution has substantial net streaming mo- 
tion and is slightly skewed in the direction opposite the 
mean velocity. 

An important idealization in the perfect ellipsoid 
model is that it asymptotes to a flat core. Most 
elliptical galaxies actually have central density cusps, 
with slopes ranging from dlnp/dr ~ 0.5 — 2.5 (e.g . 
Crane et alJ 119931: iFaber et all [19971 lLauer et all I2005L 



20071: iKormendv et all I2009D . which may be thought 



of as perturbations to the perfect ellipsoid poten- 
tial. These perturbations lead to trapping of box 
orbits by resonances (e.g. iBin nev & Tremaind I200S 
Contopoulos fe Mer tzanidca 119771: iSridhar fe Touma 



1996 : ICollett et all 119971: IMerritt fc ValimT 



Tremaine fc YuT TIoOO) , 



n ■ uj = for integer n (e.g. 
z— tube orbits are trapped by the 1:1 resonance between 
uj x and uj y ). Stars in an island of phase space around 
a stable resonance librate around the resonant orbit; 
these trapped box orbits, called "boxlets," ge n erally 
avoid the origin (e.g. iCarpintero fc Aguilarl 119981 : 
IMerritt fc Vahhirll 1199a iFulton fe Barnes! 120011 ). As 
the perturbation grows larger (the cusp gets steeper), 
more phase space is occupied by resonant islands 
until they overlap, producing regions populated by 
ergodic orbits that eventually fill the entire portion of 
their 5D energy surface not o c cupied by regular orbits 
(e.g. IMerritt fe Valluril [199& IValluri fe Merritd fl998l: 
iKandrup fe Siopisll2003l: Kalapotharakos et al.ll2004D . 

For simplicity we do not distinguish between regular 
box orbits, boxlets, and ergodic orbits in this paper, and 
will loosely use the term "box" to refer to any orbit with 
no definite sense of rotation. This will be sufficient for 
the global comparison with observations desired in this 
study. Spect ral analysis of selected orbits (Hoffman 2007; 
Clozcl 2008) reveals that the orbits classified as boxes 
are generically centrophilic (at least in the inner parts) 
and stochastic, but they remain fairly localized in phase 
space over a Hubble time. We defer a detailed study of 
the nature of the box orbits in the remnants to future 
work. 

The effect of an inner cusp on the orbital pop- 
ulations can be understood by means of a simple 



experiment - adding a central point mass to a flat-cored, 
triaxial potential with a large population of box or- 
bits. This "experiment" has been studied extensively, 
since it has direct astrophysical relev ance to black 
holes in the nuclei of g alaxies [e.g. | Norman et alJ 
19851: IGerhard fc Binneyl [19851: IQuinlan et alJ 1 1995 : 

H997; 
1998: 



Sigurdsson et al.l 119951 IMerritt fc Valluril 
Merritt fc Quintan! 119981: IValluri fc Merritti 



Poon fc Merritti 120011 : iHoll ev-Bock elmann etail 12002 



Kalapotharakos et al T l200l . Since box orbits pass 
arbitrarily close to the origin, they are deflected by the 
point mass at pericenter passage, causing the orbits to 
diffuse within box phase space. The velo city changes are 
prim arily along the axis of approach ( Chandraskchar 
1943), the x— axis for box orbits librating about the 
long axial orbit, so the angular momentum diffusion 
is mostly in the y — z plane. Eventually the star may 
wander into z— tube phase space (there is no space for 
y— tubes since they are unstable). Since box orbits also 
strike the zero-velocity surface, they convey information 
about their pericentric evolution to large radii. 

As box orbits diffuse across the z— tube boundary, 
the shape of the potential also becomes more oblate, 
shrinking the phase space availa ble for boxes and ex- 
pand ing that for z— tubes (e.g. iKalapotharakos et al.l 
l2004f h A point mass as small as 2-3% of the total mass 
can s eed this transformation in a flat-cored, triaxial sys- 
tem (IGerhard fc Binneyl 119851: IMerritt fc Valluril 119971 : 
IMerritt fc Quinlanlll998TlValluri fc Merrittlll998i r When 
the central mass concentration (CMC) is not an ideal 
point mass, the diffusion timescale varies with the degree 
of central concentration - in an r _1 cusp the timescales 
are typically longer than the lifeti me of a galaxy, 
while in an r~ 2 cusp they are short (|Merritt fc V alluri 
19961: IMerritt fc Fridmanll 1996l:lFridman fc Merrittl ll997l: 
Holley-Boc kclman n et al.l 120011 ) . In hierarchical struc- 



ture formation, there is an ongoing exchange between 
processes that induce triaxiality (e.g. violent relaxation 
in mergers) and gas inflows tha t deepen the central po- 
tential well (e.g. iDubinskil [19 94) . This process therefore 
undoubtedly plays an important role in galaxy evolution. 

1.3. Gas-rich mergers 

In current semi-analytic models, ~ 70% of the stellar 
mass in present-day ellipticals and classical bulges a s- 
sembles t hrough major mergers dHopkins et al.ll2009al P: 
iFakhouri fc Mai 120081 : iConrov fc Wechslerl 12009^ An 
elliptical galaxy-sized halo has on average undergone 
~ 1 major merger since z ~ 2 — 3, the epoch dur- 
ing which most of its stel l ar mass formed (iKauffmannl 
| 1996t iDe Lucia et al.ll2006l: iHopkins et al.ll2008al l2009al: 
iStewart et al.l P2008r i . The last major merger was typ- 
ically between two spiral galaxies, with gas fractions 
ranging from ~ 10% for systems with stellar masses 
around 3 x 10 n M w to ~ 50% for 1O 1O M systems 
(jStewart et all l2009al lErb et all 120061 and references 
therein) . 

The hypothesis that ell iptical galaxies form thr ough 
mergers between spirals (|Toomre fc Toomrel 119721 ) ac- 
tually preceded the acceptance of the concordance cos- 
mology by about two decades, based on the proper- 
ties of the galaxies themselves. Gas-rich tidal tails, 
and rings and shells indicative of the recent disrup- 
tion of a dynamically cold system, often surround 



(> 



galax ies otherw is e resembling ordin ary giant ellipti- 
cals (|Arpl 119661 : Ivan Dokkuml 120051 ). Early simula- 
tions of mergers between disk galaxies could explain a 
wide variety of the properties of observed ellipticals, 



including their r 1 / 4 l aw density pro fi les (Ivan Albada 
1981 IMcGlvnnl Il98l IHernouisti Il99l INaab fc Truiiiic 



20061 ) slow rotation and anisotropic velocity distribu- 
tions (lOstriker fc Pe ebles 19H lAarseth fc BinnevlH978 ' 



White 1978; Gerhard 1981), 



1978t iFarouki fc Shapiro! H9 82; Efstath iou et al 
fine structure (jHernquist fc Spergell 1199211 and 



flat rotation curves flWhitc 

mm, 



ap- 



parently triaxial shapes ([Gerhard! 119811 : I Barnes! 119881 : 
IFranx et allll99ll ). A simple counting argument based 
on the numbers of observed interacting pairs and ellipti- 
cal galaxies in the local universe made the prosp ect that 
these pairs turn into ellipticals quite plausible (jToomrd 
Il977l ). 

It was apparent that dissipation must play a 
large role in mergers long before hydrodynamic 
simulat ions with realistic g a s fra ct ions becam e fea- 
sible dNegroponte fc White! 119831: IHernauisti 
Nieto et alj|1991t IBarnes fc Hernquisti 119911: iLutzl 

119921 : IHernquist et all 119931 ). The 
space densities at the centers of el- 



1989; 



1991 



Ashman fc Zepl 
measured phase 
liptical galaxies far exceed the maximum densities in 
observed spirals, implying a violation of Liouville's 
theorem unless t he initial disks contain >25-30% of 
their mass in gas (!Carlberglll986t iHernquist et al.l fl993t 
iRobertson et al.ll2006af) . Dynamically cold components, 
such as embedded disks and KDCs, are often observed 
in elliptical galaxies (e.g . IFranx fc Illingworthl 119881: 



Jedrzeiewski fc Schechterl 119881: 
19911: IRixfc White! 119901: IScorza 



Hernauist & Barnes 
Bcnderl 119951: 



Jesseit et al.1 12007ft . The high specific frequencies of 



GCs in ellipticals relative to spirals imply that mergers 
must trigger the f ormation of many new clusters from 
the a vailable gas (jAshman fc Zepfl Il992l iFall fc Zhanj 
[200lh . 

More recent simulations have shown that merg- 
ers with ~ 30% gas produce remnants that fall 
on the observed fundamental plane scaling rela- 
tion (iDiorgovski fc Davis! 119871: iDressler et all 119871: 
Gonzalez-Garcia fc van Albada! 120031 : IRobertson et al.1 
2006a; Hopk ins et al.1 l2008d ). and that the remnants 
of 1:1 mergers between 40% gas disks match the ID 
kinematic properties of observed elliptical s far bet- 
ter th an dissipationless merger remnants (Co x et all 
l2006a| ). 2D kinematic maps of gas-rich merger rem- 
nants display many of the intriguing features seen in 
real galaxies, including misaligned rotation, central ve- 
locity dispersion dips, counter-rota t ing disks, and KDCs 
dJesseit et al.l 120071: IBarnes! Il992t IHernquist fc Barnes! 
11991!) The shapes of the LOSVDs of simulated 
gas-rich merger remnants display the same trends 
as el li pticals in the SAURO N sample (iBender et all 
1991 iBendo fc Barnes! l2000l INaab fc Burkertl 1 2001 1 : 
Bournaud et al.l 12005b!: iGonzalez-Garcia et all 120061 : 
Naab et all 12006a!: Hoffman et all l2009al) . and they 
occupy essentially the same part of the anisotropy- 
ellipticity (i5 — e) plane as the SAURON galaxies (with 
the notable exception of some giant slow ro tators that are 
very round, anisotropic, and featureless; iBurkert et all 
2008). These results motivate studies of the intrinsic 



orbital structure of gas-rich merger remnants, to see to 
what extent the 3D distribution functions of elliptical 
galaxies can indeed be explained with binary mergers, 
and shed light on how the observable features arise phys- 
ically. 

To our knowledge. IBarnes! (119921) was the first t o apply 
orbita l analysis in the style of Ide Zeeuwl ()1985f ) ; iStatlerl 
(|1987l ) to simulated merger remnants. He ran a series 
of dissipationless disk galaxy mergers with varying disk 
orientations and impact parameters, and classified the 
stellar orbits in the remnants based on the sign changes 
in their angular momentum. He found a wide variety 
in the remnant shapes and orbital structure, depending 
on the encounter parameters. Some of the rem nants dis- 
playe d substantial orientation twists (see also iGerhardl 
1983). The mergers often produced both x— tube and 
z— tube populations with significant net rotation, result- 
ing in large kinematic misalignments. Since the rotation 
of the majority of observed ellipticals is well-aligned with 
the major axis, he argued that dissipationless disk merg- 
e rs cannot be the gene ri c way to form elliptical galaxies. 

IBarnes fc Hernquisti (|1996f ) classified the stellar or- 
bits in mergers with 10% gas, and found that even 
this small gas component has a dramatic effect on the 
remnant shapes and orbital structure. Gravitational 
torques during the merger drain much of the gas of its 
angular momentu m, causing it to collapse inward and 
form a dense C MC (Hcrnquist 1989; Barnes fc Hernquisti 
[1991: iMihos fc Hernquisti Il994al|bl . fl996t iHopkins et all 
2009c), essentially a point mass to stars at lR e . The 
CMC destabilizes box orbits, leading to a global trans- 
formation of the remnant to a more oblate shape. 
This res ult is not surprisi n g in li ght of the studies 
by e.g . iGerhard fc Binnevl (|1985l ): iMerritt fc QuinTanl 
(1998), showing that a central point mass with just 2% 
of the total mass can globally transform the structure of 
a triaxial system. 

Uesseit et all (|2005l ): INaab et all (|2006al) performed a 
detailed study of the orbital structure of a large sample of 
equ al and unequal mass mer ger remnants with and 10% 
gas (Naab & Burkcrt 2003), with an emphasis on relat- 
ing the intrinsic structure to photometric and ID kine- 
matic observables. They studied the relation between 
the or bital structur e and i sophotal shapes in further de- 
tail in Uesseit et all (I2008D. Using spectral classification 
(|Carpintero fc Aguilarll 19981) , they found that box orbits 
typically dominate the inner parts of the dissipationless 
remnants, while x— tubes and z— tubes become dominant 
at larger radii. The box to z— tube ratio was the primary 
determinant of kinematic properties such as the location 
of the remnants in the h^^ — v/a planes. When a gas 
component was added, the box population was highly 
suppressed and the remnants became z— tube dominated 
and oblate. The shape of the z— tube orbits in the more 
axisymmetric dissipative remnants made the isophotes 
less boxy. The 1:1 merger remnants were slowly rotating, 
while the 3:1 remnants were found to be rapidly rotating 
and disky. They concluded that observed rapidly rotat- 
ing ellipticals could form from di ssipative 3:1 mergers, 
but b oxy, slowly rotating systems (jKormendv fc Bender! 
1996) could not have formed through dissipative disk 
mergers. 

In this work we analyze the orbital structure of 1:1 
merger remnants in simulations including SF and feed- 
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back, enabling us to consider the hi gh gas fractions char- 
acteristic of spiral galaxies at z ~ 2 ([Stewart et al.ll2009at 
iDaddi et al.l 120091 ). At fixed f gas the effect on the rem- 
nant structure may be diminished by including SF and 
the dissipational features may become more spatially ex- 
tended, since the gas can be converted to collisionless 
material early on in the merger and undergo subsequent 
violent relaxation. We quantify the variation of the or- 
bital structure and intrinsic shape with gas fraction for 
fgas ranging from to 40%, and discuss the physical 
mechanisms that may be driving the orbital transforma- 
tions, with direct comparisons against dynamical mod- 
els of observed systems in mind. We relate the orbital 
structure of the remnants to their appearance in 2D kine- 
matic maps on lR e and 3i? e scales, and show that a wide 
range of different kinematic structures can be accounted 
for simply by varying f gas in 1:1 disk mergers. 

1.4. Outline 

In §2 we describe our merger simulations and remnant 
analysis methods, and in §3 we present our results. The 
intrinsic structure of our eight dissipationless remnants 
is presented in §3.1, as a baseline for understanding the 
effect of adding a gas component. We show how the or- 
bital structure and shape of the remnants varies with f gas 
in §3.2, and discuss the KDCs and embedded disks that 
arise in the gas-rich remnants in §3.3. In §3.4 we briefly 
describe the structure of the remnants beyond ~ li? e , de- 
ferring a more detailed and quantitative presentation to 
iHoffman et ail (|2009bD . In §4 we summarize our results 
and conclude. The radial profiles of the orbital struc- 
ture, intrinsic shape, and orientation of all 48 dissipative 
remnants are presented in the Appendix. 

2. SIMULATIONS AND METHODS 

2.1. Simulations 

Our galaxy merger simulations were performed us- 
ing the publicly available TreeSPH ( Smoothed Par- 
ticle Hydrodynamic s) code Gadget-2 ([Springel! 120051 : 
ISpringel et all 120011) . which uses an advanced formula- 
tion of SPH that explicitl y conserves both energy an d 
entropy when appropriate (Springcl fc Hernquistl 12002). 
In addition to the standard features, our version of the 
software inc ludes sub-resolut i on prescriptions for radia- 
tive cooling (|Katz et al.lll996HDave et al.lll999ft . SF, and 
feedb ack from supernovae an d AGN, as described in de- 
tail in ISpringel et al.l (|2005a|) . 



SF and supernova fe edback are treated using the m ulti- 
phase prescrip t ion of ISpringel &; Hernquistl ((20031 ) and 
ISpringel etakl (|2005aD . The interstellar medium (ISM) 
consists of a cold cloud phase in which stars form, and 
a hot phase that provides pressure support for the disk. 
The density dependence of the SF rate is cali brated to 
matc h the observed S chmidt-Kennicutt Law (| Schmidt! 
119591 : iKennicuttl [1998). Radiative cooling drives gas 
transfer from the hot phase to the cold phase, while heat- 
ing from supernovae feeds gas back into the hot phase. 
The parameter qeos smoothly varies the effective equa- 
tion of state (EOS) between isothermal {qEOS — 0) and 
the pure multi-phase model (qEOS = 1); higher values 
of qEOS correspond to a stiffer EOS and permit stable 
disks at higher gas fractions. In all of the simulations in 
this paper we used qEOS = 0.25. 

The simulations also include sink particles repre- 
senting black holes, which accrete gas at a rate 

M B h = mm(M B ondi, M E dd), where M Bon di oc 
M BH p gas / yjc 2 s + v 2 is a rate bas e d on the Bondi-Hoyle- 
Lyttleton formula §Bon M SHj iBondi fc Hovli IT941 
iHovle fc Lvttletonl [19391) . and M E dd is the Eddington 
rate. A fraction tft r of the accretion energy is re- 
leased thermally and isotropically into the surrounding 
gas, where e r = 0.1 is the assumed radiative efficiency, 
and ef — 0.05 is the fraction of the radiated luminosity 
that can couple thermally to the gas. 

The gravitational softening length was e = 140 pc for 
the stellar and gas components in our simulations, and 
e = 290 pc for the halo component. We decreased the 
parameter rj controlling the Gadget-2 timestep criterion 
(At = ^/2r]e/\a\, where a is the particle acceleration) 
from the default value of 0.025 to 0.0025, to ensure that 
the steep central cusp that formed in the gas-rich simu- 
lations was preserved for several Gyrs after the merger. 

We set up our merger ICs u s ing st andard tech- 
niques presented in ISpringel et al.l (|2005aD . which have 
been emplo yed in a varie ty of previous applica- 
tions (e.g. iDi Matteo et al.l 120051: iNaab et all [2006 
Hopkins et al.l 120061: iRobertson et al.l l2006at ICox et al 



2006aH Johansson et al.ll2009HHopkins et al.l l2009b). The 
galaxy models consisted of exponential disks embedded 
in dark matter halos; stellar bulges were not included 
in this study. The total mass was specified through the 
virial velocity, ^200 = 160 km/s for an approximately 
Milky Way-sized gala xy, by M tot = vL n /(10GH Q ). The 
halo was modeled as a lHernquistl ([ 19901 ) profile with spin 
parameter A = 0.033, and scale radius chosen to match 
an NFW profile with concentration c = 0.9 in the inner 
parts. 

The disk mass was set to Md — fdM to t, with fd — 
0.041. Its scale length of = 3.9 kpc was determined by 
the requirement that it be centrifugally supported with 
the same specific angular momentum as the halo, Jd — 
fdJtot- The stellar component of the disk was assigned a 
radially constant vertical scale height of 0.2rd, while the 
vertical structure of the gas component was set by the 
requirement of hydrostatic equilibrium. The disk models 
were realized with 80000 equal-mass particles, a fraction 
1 — f gas of them in collisionless stars and the other f gas 
in SPH particles. The halo was comprised of 120000 
collisionless dark matter particles. 

For each simulation two identical disk galaxies were 



constructed following this procedure, and placed on a 
parabolic orbit with an impact parameter of 7.1 kpc. At 
their initial separation of 140 kpc, the center-of-mass tra- 
jectories were well- approximated by point mass orbits. 
At least in the inner parts, the final structure is not 
highly sensitive to our choice of small impact parame- 
ter, since the most bound particles (the stars) lose most 
of their orbital angular momentum to the halo by the 
time the cores merge even in wide encounters ([Barnes! 
1998: ICox et afl l2006a) . The effect of the merger impact 
parameter on t he stellar halo structu re will be explored 
in future work ([Hoffman et al.H2009bfl . 

The orientation of the disks relative to the merger 
(x— y) plane was parametrized by the two angles 9 and <f>, 
the polar angle of the spin vector relative to the z— axis 
and its azimuthal angle relative to the axis o f app roach, 
as shown in Figure 6 of iToomre fc Toomrd (|1972t ). We 
used the scheme proposed bv lBarnei (|1992f ) to sample the 
space of possible disk orientations in a relatively unbiased 
way, including prograde, retrograde, and polar orbits. 
The spin of the first disk coincides with each of the four 
symmetry axes of a regular tetrahedron pointing upward, 
while that of the second disk points along the axes of the 
corresponding downward-pointing tetrahedron. The re- 
sulting eight merger orbits are listed in Table 1. This set 
of encounter orbits has been used in a number of pre- 
vious studi es (e.g lBarnesHl992t iNaab fc Burkertl [20031: 
INaab et alj|2006at ICox et al.H2006al) . which may be used 
as a reference point for comparison with our results where 
appropriate. The eight encounters listed in Table 1 were 
repeated at seven different gas fractions, f gas = 0, 5, 10, 
15, 20, 30, and 40%, for a total of 56 merger simulations. 

Once on their orbits, the galaxies reach their first peri- 
center passage at i « 0.35 Gyrs, at which point their 
morphologies become strongly distorted by tidal forces 
from the other galaxy, and the stars and gas each form a 
bar. The stellar bar slightly trails the gas bar, torquing 
back on the gas and draining its angular momentum 
so that the gas flows inward and undergoes a burst of 
star S F (iMihos fc Hernauistlll994at iBames fc Hernquistl 
fl¥96t lEscala et all 12004 iHopkins et alJl2009b[ ). At t » 
1.8 Gyrs the galaxies reach second passage, and their 
cores merge shortly thereafter. The strong and persis- 
tent gravitational torques during the second passage and 
final merger typically induce a stronger central starburst 
than at first passage. By about 0.5 Gyrs after the merger 
the size, shape, and velocity dispersion of th e remnant 
(meas ured within ~ lR e ) reach a steady state ([Cox et al.l 
2006a|), and the system may be considered relaxed. 

2.2. Remnant analysis 

We freeze the potential at t = 4.3 Gyrs (about 
2.5 Gyrs after the merger is complete), and represent 
it using the bi-orthogonal "designer" basis expansion 
(IClutton-Brockl U97l Il97l IBinnev fc Tremainel [20081) 
presented in iHernquist fc Ostrikerl (|1992f ). 

p(r, 9, 4>) = ^ A nlmPnl(r)Y lm (e, 4>) (1) 

nlm. 

$(r, 9, 4>)=J2 B mm$ni(r)Yi m (0, <t>). (2) 

nlm 

The terms in this expansion individuall y satisfy Po i sson's 
equation, and the lowest-order term is a lHernquisti (1990) 



profile, which should be a good first approximation to 
the density profile of a remnant resembling an elliptical 
galaxy. Before computing the expansion coefficients we 
symmetrize the particle distribution by reflecting all of 
the positions and velocities about the origin, effectively 
increasing the particle statistics by a factor of two and 
reducing global fluctuations. The orbit of each stellar 
particle is then integrated through ~ 300 radial turning 
points (peri/apocenter passages) in the static potential 
using a Bulirsch- Stoer integrator ([Bulirsch fc Stoei1ll966t 
IPress et al.lll992l) . This technique allows all of the stellar 
orbits to be efficiently followed for ^150 dynamical times, 
including halo orbits for which this corresponds to many 
Hubble times. It also lessens spurious relaxation and 
simplifies the orbital analysis by fixing the principle axes 
of the system. 

Each stellar orbit was classified as a box, a:— tube, or 
z— tube orbit using a simple algorithm based on the sign 
changes in th e star's angular momentum, introduced by 
I Barnes! (|1992t ). At each timestep we computed the angu- 
lar momentum vector j, and checked whether each com- 
ponent of j had changed sign since the last step. At 
the end of the integration we constructed the vector k 
such that hi = 1 if ji never changed sign, and ki — 
otherwise. The orbit was assigned the classification code 
C = k x + 2k y + Ak z , which is for a box orbit, 1 for an 
a:— tube orbit, and 4 for a z— tube orbit. For a perfect 
triaxial ellipsoid, other values of C are unphysical. 

This classification scheme depends on a correct speci- 
fication of the system orientation, and many of the rem- 
nants have intrinsic orientation twists. We therefore di- 
agon alized the inertia tensor, in the form 7y = XiXj 
(e.g. lAguilar fc Merrittl 119901) . in ~ 50 cumulative en- 
ergy bins, and used the local orientation based on the 
star's energy for the classification. The eigenvectors of 
Iij are the principal axes of the remnant's figure, and 
the square roots of the eigenvalues give the relative scale 
lengths along the three principal axes, a, b, and c. All 
particles - gas, stars, and dark matter - were included in 
Iij, since a star's orbit is determined by the sum of all 
three components. Note that since we used the local ori- 
entation for the classification, e.g. the position angle of 
the z— tube rotation may appear to vary with galactocen- 
tric radius in the remnants with significant orientation 
twists. 

Only ~1% of the stellar orbits were not assigned C 
values of 0, 1, or 4, except in rare instances of sud- 
den orientation twists or locations where the potential 
was very nearly spherical. The results of our simple 
classificat ion scheme were consistent with spectral clas- 
sifica tion (iCarpintero fc Aguilarll 19981 : IBinnev fc Sperge] 
Il982t 11984 lHoffmanll2007l ). and we determined that the 
simpler algorithm was better suited for the noisy poten- 
tials and gross analysis desired in this work. However 
note that this algorithm does not distinguish between 
boxes, resonant boxlets, and ergodic orbits. 

For each remnant, we present radial profiles of the or- 
bital structure and intrinsic shape as shown for one ex- 
ample in Figure [3l First we order the stellar particles 
by energy, and divide them into ~ 20 energy bins with 
equal numbers of particles. The mass fraction of stars 
in box, x— tube, and z— tube orbits is computed for each 
energy bin, and plotted vs. percentile in binding energy 
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Fig. 3. — Radial profiles of the orbital structure and intrinsic shape of remnant i, 40% gas. Left: Fraction of the stellar mass 
in each of the three major orbital classes vs. percentile in binding energy. The horizontal axis is labeled by the mean radius (in units of 
R e ) of the stars in the current energy bin (see text for details). Left center: Rotation bias, p = (JV+ — N—)/(N+ + N-), of the a;— tube 
and z— tube orbits, ft = corresponds to perfect cancellation of the rotation, while j3 = ±1 means that all of the stars are streaming in 
the same direction. We choose the overall sign so that /3 is positive at the 40th percentile in binding energy (lR e ). Right center: Intrinsic 
shapes of the remnants. The heavy green line is the triaxiality parameter, T, for all particles in the current energy bin and more bound. 
The black line is the intrinsic ellipticity e = c/a, and the red dashed line is the smaller of the other two axis ratios. The horizontal line 
at T = 0.5 marks the boundary between oblate and prolate shapes. Both stars and DM particles are included in the evaluation of the 
intrinsic shapes. Right: Intrinsic orientations (9 symm , 4> symm ) and kinematic misalignments (^kin), in degrees. The red solid and dashed 
lines are the polar (9 syrnm ) and azimuthal (tfisymm) angles of the symmetry axis, in this case the short axis since the remnant is oblate. In 
general we will plot 6, rf> of the short axis in red if the remnant is oblate (T < 0.5) at lR e , and 8, <fi of the long axis in blue if the remnant 
is prolate (T > 0.5) at lR e . 9 and are evaluated in a coordinate system where the disks' initial orbital angular momentum points in 
the 2 direction, and the galaxies initially approach one another along the x-axis. To evaluate ^kin (equation 3; black line), we return to 
the local principal axes coordinate system. In all panels, the lR e scale probed by SAURON and the 3R e scale typical of SKiMS data are 
marked with vertical lines. 



(by mass), as in the left-hand panel of the figure. For 
ease of comparison with observational results, we also 
compute the mean galactocentric radius of the particles 
in each energy bin, and label the bins corresponding to 
0.3, 1, 3, and the largest integer number of R e . 

For understanding kinematic maps, it is useful to quan- 
tify not only the populations of the major orbital classes, 
but also their degree of streaming. For this purpose we 
plot profiles of 

0= (N+-N-)/(N+ + N-) (3) 

for the a;— tube and z— tube orbits, as shown in the left 
center panel. Here N+ and AL. are the numbers of par- 
ticles rotating in either direction; /3 = corresponds to 
perfect cancellation of the rotation, while (3 — ±1 means 
that all of the stars are streaming in the same direction. 
The sign convention was chosen by requiring (3 to be pos- 
itive at the 40th percentile in binding energy (lR e ). 

The intrinsic shape and orientation vector are com- 
puted by diagonalizing the inertia tensor as described 
previously. We plot the maximum and minimum axis 
ratios, e — c/a and r m i n = mm(b/a,c/b), and the triaxi- 
ality parameter, 

T=(a 2 -b 2 )/(a 2 -c 2 ), (4) 

in the right center panel of Figure [3J T ranges from for 
a perfectly oblate spheroid to 1 for a prolate spheroid. In 
the right-hand panel we plot the polar and azimuthal an- 
gles (Osymm and 4> S ymm) of the symmetry axis - the long 
axis if T > 0.5 within lR e , or the short axis if T < 0.5. 
If the short axis orientation is plotted then we color the 
lines red; otherwise we color them blue. The shape of a 
remnant that is oblate within lR e may become prolate 
farther out; in this case we still plot the orientation of 
the short axis at all radii, which does not make the an- 



gles ill-defined since the remnants are never too close to 
axisymmetric at large radii. Note that 9 symm and 4> S ymm 
represent angles between a vector and an axis (whose di- 
rection is defined only up to a sign) , and therefore range 
from to 90°. When plotting the orientation profiles, we 
use a fixed coordinate system where the z— axis points 
along the disks' initial orbital angular momentum vector, 
and the x— axis is along the initial direction of approach. 

On the same axes we plot the intrinsic kinematic mis- 
alignment, 

*fein = arctan \j z /j x | , (5) 

or the angle between the short axis and the net angular 
momentum vector. To compute ^kin we revert back to 
the local principal axes system, so that ^kin = always 
corresponds to rotation about the intrinsic short axis. 

To place our work in an observational context, it is 
also useful to show how the intrinsic orbital structure 
imprints itself on SAURON-like 2D kinematic maps. We 
constructed histograms of the LOS velocity in 40 x 40 
spatial bins within lR e , using 80 velocity bins within 
±3.5ct. For each LOSVD we performed a least-squares 
fit to the 5-parameter function 

F(y) = Ae- w2 / 2 [l + h 3 H 3 (w) + h 4 H 4 (w)}, (6) 

where w = (y — v)/a, and v and a represent the mean 
velocity and velocity dispersion. Here H 3 and H4 are 
Gauss-Hcrmitc (GH) polynomials, and the parameters 
/13 and h 4 measure the skewness and ku r tosis of the dis- 
tribution (Ivan der Marel fc Franx! U995 IGerhardl fl993; 
IBinnev fc Merrifieldlll998Tr The Gadget particles were 
smoothed over a radius h smoot h = max(/i see , h n gb), where 
h see — 150 pc corresponds to a seeing of 1.5" at 20Mpc, 
and h ng b is 1.7 times the distance to the 128th nearest 
neighbor. 



10 



We focus primarily on the velocity fields in this work. 
Note that v, as derived from the fit of equation 4, deviates 
systematically from the true mean of the distribution for 
nonzero J13 and hi (see Ivan der Marel fc Frarudll993l for 
the relevant correction terms). Ho wever we follow the 
convention of previous authors (e.g. iBender et al.l[l994l ) 
and simply use v from the GH fit to represent the velocity 
in this paper. 

When displaying all of the remnants at a given f gas , we 
choose the viewing angles randomly (from an isotropic 
distribution) to give a fair representation of their ap- 
pearance in kinematic maps. We label each map with 
the LOS angles and <j>, where 8 is the polar angle be- 
tween the LOS and the z— axis, and <j) is the azimuthal 
angle between the LOS and the y— axis. The apparent 
major-axis rotation is maximized when 9 = 90°, and the 
minor-axis rotation is maximal at <f> = for a given 6. 

3. RESULTS AND DISCUSSION 
3.1. The dissipationless remnants 

The structure of the dissipationless remnants is pre- 
sented in Figure 01 There is substantial variation in 
the orbital distribution and shape over the eight merger 
orbits - for instance remnant j is uniformly prolate 
and dominated by x— tube orbits, while remnant o is 
oblate-triaxial and composed primarily of box orbits and 
z— tubes. Overall the remnants are dominated by box or- 
bits in their inner parts (r < 1.5R e ), and tub e orbits in 
their o utskirts. This is not surprising since the lHernquist] 
( 1990) profile, which generically provides a good fit to the 
density profiles of dissipationless merger remnants, has a 
shallow inner cusp (p oc r" 1 ) and much steeper outer 
profile (p cx r~ 4 ) Q. 

The inner shapes of the remnants are typically prolate- 
triaxial and highly flattened (T ~ 0.8, e ~ 0.6), consis- 
tent with previous results for head-on mergers of stellar 
systems that have sh ed most of their orbita l angular mo- 
mentum to t he halo dMiller fc Smitbl ll98Ct IWhitelll983t 
iBarnesI 1199$ iCox et al.l l2006al ). Farther out the rem- 
nants become roun der (e ~ 0. 8 ) and cl oser to maximally 
triaxial. See also iDubinskil (119941): iKazantzidis et alj 
(12001; iDebattista et alJ (1200811: IValluri et all (12009) ; 
lAbadi et al.l (j2009D ; IRomano-Diaz et al.l (|2009T) . who 
studied halo shapes in a cosmological context and found 
that the response of the DM to baryonic condensation 
into galaxies drives the halos to rounder and more oblate 
shapes. 

There is substantial cancellation of the prograde and 
retrograde rotation of the z— tube orbits (/3 near zero), 
especially in the inner parts. On the other hand the 
a:— tube populations are highly streaming, with j3 typ- 
ically approaching one at large radii. Remnant j is a 
notable exception to this rule, with the rotation of its 
dominant x— tube orbits canceling nearly perfectly over 
the full range of radii plotted. Upon closer examination 
of this remnant, it turns out that this owes to the sym- 
metry of the ICs. Both disks on orbit j are substantially 
inclined, but their spins point in opposite directions. If 

A star in what we are calling the "outer parts" (~ 3 — 6_R e ) still 
occupies the inner DM halo (which also follows roughly a Hernquist 
1990 profile but with a much larger scale radius). The total density 
(stars + DM) at these radii scales roughly like r~ 2 - still steep 
compared to r -1 , but not nearly as steep as r — 4 . 



the x— tube orbits in the remnant are separated by their 
disk of origin, nearly all of those from disk 1 rotate in 
one sense about the long axis, while those from disk 2 
rotate in the opposite sense, giving the net cancellation 
seen in Figure 2] 

The dissipationless remnants do not show strong ori- 
entation twists, and their long axis (which is typically 
the axis of symmetry, since the remnants are prolate) 
lies in the merger plane {9 symm — 90°), along the axis of 
the final coalescence of the stellar components. This i s 
consistent with the interpretation of lNovak et al.l {2006) : 
the time- varying tidal field along the merger axis cou- 
ples coherently to the stellar orbits, producing a bounce 
that elongate s the remnant along the direction of tidal 
compression (|Millerl 119781: IWhitel(l983l) . The kinematic 
misalignments are large since the net rotation tends to 
favor the long axis, and noisy since the magnitude of the 
rotation is often small. 

In the li? e kinematic maps, the dissipationless rem- 
nants are generally slowly rotating (e.g. remnants j and 
p) or are dominated by minor-axis rotation (e.g. rem- 
nants m and n). Some of the remnants (e.g. p and 
I) definitely show both major and minor axis rotation, 
though small in magnitude (v/a ~ 0.15). 

Figure U is a baseline for understanding how a gas 
component affects the remnants - it may be compared 
with Figures [15] - [20] in the Appendix to see how the or- 
bital structure evolves for each merger orbit as f gas is 
increased from to 40%. 

3.2. Variation with gas fraction 

In Figure [5] we illustrate the effect of dissipation on 
the remnant structure by varying the gas fraction for a 
fixed merger orbit, i, whose structural transformations 
are representative of the overall trends. Case i is an 
encounter between a prograde disk in the merger plane 
(#1 — 4>i = 0), and a highly inclined disk that is also 
slightly prograde (0 2 = 71°, <p 2 = 30°). The orbital 
structure, shape, and orientation profile of this remnant 
are shown for seven different gas fractions in Figure [5J 
Unlike in Figure UJ the velocity maps are all in projection 
along the y— axis, to clearly show the rotation of both the 
x— tube and z— tube orbits and how it varies with f gas - 

As fg as increases from to 40%, the population of box 
orbits within ~ 1.5i? e declines, and the z— tube popu- 
lation increases. The most rapid change in the orbital 
populations occurs between and 20% gas. The outer 
orbital structure (r > 2i? e ) is relatively unaffected by 
the gas, and the a;— tube population does not display an 
obvious trend. 

The variation in the intrinsic shape closely follows the 
change in the orbital populations. As f gas increases, the 
inner parts become progressively more oblate; the 30- 
40% gas remnants are nearly oblate axisymmetric within 
li? e (T < 0.1). The remnants also become rounder with 
higher f gas '- e increases from ^0.6 for the dissipationless 
remnant to ~0.8 at the highest gas fractions. 

Though the mass in z— tube orbits within ~ 2R e rises 
substantially as f gas increases from to 15%, the rem- 
nants show little change in their major-axis rotation over 
this range, since the mass in prograde and retrograde or- 
bits is nearly equal (f3 ~ 0). The 15% gas remnant, 
though dominated by z— tubes within lR e , still appears 
as a minor-axis rotator in the projected velocity map. 
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Fig. 4. — Intrinsic structure of the dissipationless remnants. The first four columns follow the format of Figure [3] The merger 
orbit (see Table 1) is specified on the left for each row. In the right-hand column we show a lR e velocity map of the remnant, viewed along 
a randomly (isotropically) chosen LOS. The two angles specifying the LOS are indicated by the labels on the right. 8 is the polar angle 
between the LOS and the short axis (this axis was chosen as the reference throughout the paper since most of the remnants with gas are 
oblate), and <f> is the azimuthal angle between the LOS and the intermediate axis. The major-axis rotation is maximal at 8 = 90°, and the 
minor-axis rotation is maximized at tj> = for fixed 8. The velocity scale runs from — v max to v ma x, where v max is the maximum value of 
\v\. The value of v max (in km/s) is indicated by the label in the upper right corner of each map. 
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Fig. 5. — Evolution of the intrinsic structure of remnant i with gas fraction. The format is the same as in Figure [4] except 
that now each row corresponds to a different f ga s- In this figure all of the velocity maps are shown in projection along the intermediate 
axis, (9, <j>) = (0, 0), to allow the major- and minor-axis rotation in the remnants at different f ga s to be compared directly. 



13 



v (km / s) 



-70 70 




-450 -150 150 450 -450 -150 150 450 



v (km / s) v (km / s) 

Fig. 6. — Cancellation of the z— tubes and streaming of 
the x— tube orbits. Upper left: Velocity map of remnant m, 10% 
gas, viewed in a near-intermediate axis projection. Upper right: 
Same map with x— tube orbits removed. Lower left: LOSVD at 
location A, near the major axis. Lower right: LOSVD at location 
B, near the minor axis. The velocity scale is the same on both 
maps. 

In the 20% gas remnant, the innermost z— tube orbits 
begin to show substantial streaming, in a sense opposite 
the net rotation farther out. This inner streaming shows 
up as a prominent KDC in the velocity map. Note that 
the minor-axis rotation present at lower gas fractions is 
still there in this map, but it is dwarfed on the veloc- 
ity scale by the rapid rotation within the KDC. In the 
30-40% gas remnants, the z— tube orbits show significant 
streaming (f3 ~ 0.3) over the full range of radii plotted, 
and the edge-on kinematic maps display rapid, disk-like 
major-axis rotation. This rapid rotation about the in- 
trinsic short axis gives kinematic misalignments near zero 
in the inner parts. 

Although the rotation bias of the z— tube orbits in the 
30% gas remnant appears fairly uniform, the cause of the 
streaming in the inner and outer parts is quite different. 
In the outer parts the net angular momentum is retained 
from the ICs, and is therefore present at all gas fractions. 
The inner streaming is a direct consequence of dissipa- 
tion - the gas that retains its angular momentum during 
the merger re-forms a cold disk as it dissipates energy, 
which forms new stars primarily after the violent relax- 
ation of the merger is complete. Thus the inner part of 
the gas-rich remnants is a superposition of a hot, slowly- 
rotating component consisting of the stars present before 
the merger (hereafter referred to as "old stars"), and a 
maximally rotating disk composed of "new stars" formed 
through dissipation. This is why /3 is only around 0.3 de- 
spite the high maximum velocities in the maps. We will 
elaborate on this point in section 3.4. 

The inner parts of the gas-rich remnants are often 
rapidly tumbling, which produces large intrinsic orien- 



tation twists between ~1 and 3i? e , as in the bottom two 
rows of in Figure [5l Note that the orientations in the 
inner parts of these remnants are an arbitrary snapshot 
in time; had the simulation been frozen earlier or later 
Qsymrn and (ftsymm could have quite different values. 

In the preceding discussion we have demonstrated that 
the low- / sa ,5 remnants show little major-axis rotation not 
because they lack z— tube orbits, but because there is a 
large degree of cancellation between the z— tubes rotating 
in opposite directions. An observer who measured the 
major- and minor-axis rotation speeds of the typical 10% 
gas remnant might naively think that the system was 
prolate in shape and dominated by a;— tube orbits, when 
in fact the predominance of the minor-axis rotation is 
caused by a sub-dominant population of x— tubes that 
are highly streaming. 

We highlight this point in Figure |6] showing remnant 
m at 10% gas, which has comparable masses in x— tubes 
and z— tubes within li? e , but somewhat more z— tubes 
(see Figure [T6f . The upper left panel is a velocity map 
of the remnant in a nearly edge-on projection (maximiz- 
ing both the major- and minor-axis rotation), and the 
second panel shows the same map with the x— tube or- 
bits removed. Without the a;— tubes, the remnant shows 
almost no rotation at all. The LOSVDs at the locations 
marked by the white letters are plotted in the lower two 
panels. The lower left panel clearly shows the canceling 
streams of z— tube orbits along the major axis, while the 
lower right panel depicts the highly streaming x— tubes 
along the minor axis, peaked at nearly 200 km/s. In the 
lower left panel, the large population of box orbits fills 
in the gap in the bimodal z— tube distribution, yielding 
a total LOSVD that is nearly Gaussian in shape. 

Figure [7] provides insight into the physical explanation 
for this orbital structure by mapping out the locations, in 
the original disks, of the stars that end up in each of the 
three orbital classes in remnant i, with f gas =0, 15, and 
30%. Disk 1 is the direct disk, and disk 2 is the inclined 
one. Only old stars are included in this plot, so the dissi- 
pative disk component in the 30% gas remnant is left out. 
Each disk was partitioned into cylindrical bins contain- 
ing 300 stars apiece, and the bins were labeled accord- 
ing to which orbital class formed the plurality of stars 
among those 300 in the final remnant - black for boxes, 
green for z— tubes, and red for cc— tubes. The bins fa- 
voring tube orbits were further labeled based on whether 
they tended to stream in one direction in the final rem- 
nant, or to have canceling rotation. Bins with /3 > 0.3 
are indicated with upward-pointing arrows and labeled 
as "counterclockwise," while those with f3 < —0.3 are 
indicated with downward-pointing arrows and labeled as 
"clockwise." Bins with \f3\ < 0.3 are labeled as "cancel- 
ing" and indicated with filled circles. The blue circles 
mark the half-mass radii of the original disks. 

In the dissipationless simulation, most of the stars 
within the half-mass radius become box orbits in the fi- 
nal remnant. When a dissipative component is added 
(second and third columns), the set of stars that be- 
came box orbits in the collisionless case instead form 
a population of z— tube orbits with largely canceling 
rotation, especially in disk 2. It is apparent that 
this cancellation does not arise primarily from the ICs, 
since it occurs for stars in the same disk. Note that 
the spatial boundary of the canceling z— tube popu- 
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Fig. 7. — Spatial distribution in the original disks of stars that end up on box, a;— tube, and z— tube orbits, for merger 
orbit i. Disk 1 (top row) is on a direct, prograde orbit (9\ = <t>\ = 0), while disk 2 (bottom row) is inclined and prograde ($2 = 71°, 
<p2 = 30°). Each disk was partitioned into cylindrical bins containing 300 stars each, and we identified the dominant orbital class of the 
stars from each bin in the final remnant. If tubes were dominant, then we also computed /3 for the bin. Bins with \/3\ < 0.3 are labeled 
with circles, while those with |/3| > 0.3 are labeled with triangles pointing up or down according to the sense of rotation (clockwise or 
counterclockwise, where the overall sign is arbitrary). The key for different orbital classes is given in the legend. The three columns show 
the remnants at 0, 15, and 30% gas. Only old stars are included (the gas particles in the initial disks are excluded). The blue circles mark 
the half-mass radii of the original disks. 



lation changes little between 15 and 30% gas, sug- 
gesting that the boundary is set by a bar forma- 
tion criterion that determines which orbits "would have 
been" boxes (e.g. iMiller fc Smith! [1971 iNorman fc Silkl 
19831: iBarnes fc Hernauistl 119961: iBournaud et al.l [2005a; 
MacMillan et all 120061: IHopkins et al.l l2009bTT ^ rather 
than by the direct gravity of the dissipative compo- 
nent. The box population might be fully transformed 
into z— tubes so long as a critical gas mass is exceeded. 
It should also be noted that the efficiency with which the 
gas is torqued inward decreases with increasing f gas , so 
the CMC mass comprises a smaller fraction of the ini- 
tial gas mass in high er- f oas remnants (see Figure 7 of 
IHopkins et afll2009bl) . 

Stars outside the half-mass radius in disk 1 generally 
become z— tube orbits in the remnant, with a higher de- 
gree of streaming than in the inner parts. The fate of 
these stars is relatively insensitive to f gas , though some 
box orbits from this region in the dissipationless rem- 
nant do become z— tubes in the simulations with gas 
(compare the top right and top middle panels). The 
stars originating in the outer part of disk 2 end up in a 
uniformly streaming x— tube population, irrespective of 
fgas- Their streaming about the long axis of the remnant 
can be explained simply in terms of the mer ger geome- 
try - see e.g. Figure 3 of iNovak et all (f2006h . The long- 
axis lies in the merger plane, so stars on nearly circular 



orbits in the outskirts of an inclined disk start out with 
large tangential velocity components about this axis, and 
are t herefore liable to occupy th e ie— tube "start space" 
(e.g. lSchwarzschiidlll979L Il99l . The uniform rotation 
of these orbits in the final remnant suggests that they do 
not cross orbital boundaries during the merger; their an- 
gular momentum is retained from the initial spin of the 
inclined disk, and is a direct signature of the dynamically 
cold nature of the progenitor galaxy. Figure [7] provides 
a good example of how orbital classification can isolate 
sets of stars with similar physical histories. 

So far we have focused the discussion on individual 
remnants; Figure [5] shows the dependence of the orbital 
structure on f gas , averaged over all eight merger orbits. 
We compute the structural parameters on ~ lR e scales 
since this radius has been the main fo cus of triaxial dy- 
namical modeling efforts to date (e.g. Ivan de Ven et al.l 
2008t Ivan den Bosch et aD 120081: Ide Lorenzi et al.l l2006t 
van den Boschll2008f ) . and defer further discussion of the 
outer orbital structure to section 3.5. 

As expected from the preceding discussion, the z— tube 
population increases, the mass in box orbits decreases, 
and the shapes become rounder and more oblate with 
increasing f gas - This transformation occurs primarily 
between and 15% gas. At higher gas fractions the av- 
erage shapes, box, and z— tube fractions level off, con- 
sistent with the "critical mass" interpretation of the box 
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Fig. 8. — Inner structure (< IRe), averaged over eight merger orbits. Upper left: Fraction of the stellar mass in each orbital class 
vs. fgas, for the 25-45th percentile in binding energy (~ 0.5 — 1.3i? e ). Upper right: Rotation bias (/3) of the x— tube and z— tube orbits, as 
denned in equation 1. Lower left: Intrinsic shape of the 40% most bound particles. The filled diamonds are for all matter (including DM 
particles), while the asterisks are for stars only. Lower right: Intrinsic orientation twist between the 25th and 50th percentile in binding 
energy, and kinematic misalignment (^kin) f° r the 40% most bound particles, ^kin is defined as in Figure [3] The orientation twist is 
computed as max(Bij), «25 < (*>j) < *50> the maximum angle between the symmetry axes computed within any two cumulative energy 
bins between the 25th and 50th percentiles. 



orbit suppression. The apparent break in the trend at 
30% gas arises because these remnants are very round. 
In a remnant that is very nearly spherical, the distinction 
between the two classes of tubes is blurred, and the tri- 
axiality parameter and orientation vector are ill-defined 
(see e.g. remnant m in Figure H"9|). The 40% gas rem- 
nants are slightly more flattened than those at 30% gas, 
presumably owing to the rotation of the stronger em- 
bedded disks. The stellar component is slightly more 
flattened than the mass including the DM, particularly 
at low f gas . There is a mild downward trend in the pop- 
ulation of x— tube orbits with f gas between zero and 15% 
as the shape of the potential becomes more oblate, but 
the round potentials at higher gas fractions are slightly 
more amenable to a;— tubes. 

At low fg as the x— tube orbits tend to stream far more 
than the z— tube orbits. This trend reverses at the high- 
est gas fractions owing to the rapid rotation of the em- 
bedded disks, and the direct modification of the gravita- 
tional potential by the large CMC, which may blur the 
lines between x— tube and z— tube orbits as the potential 
becomes very round. The streaming of the a;— tube or- 



bits might also fall off because of the figure rotation; the 
long axis no longer necessarily lies in the merger plane, 
so more stars may be randomly scattered onto tangen- 
tial orbits about the long axis. The exact reason for 
this trend in the x— tube rotation needs further investi- 
gation. The cancellation of the z— tube rotation is most 
complete around 10% gas, perhaps because the stronger 
CMC more efficiently disrupts box orbits than at lower 
fgas, while at higher f gas the dissipation begins to play 
a direct role. 

In the high-/ gas remnants, the fraction of box orbits 
often turns upward in the very center. This effect is 
partially hidden in e.g. Figure [20] since we started the 
profiles at the 15th percentile in binding energy. It can 
be naturally explained since only the innermost compo- 
nent (the CMC itself) contains no smaller CMC, and so 
can support stable box orbits. The violent relaxation 
during the final merger of the two cores provides a nat- 
ural setting for the onset of triaxiality. However the box 
population computed within the central starburst is ex- 
pected to be highly sensitive to numerical effects, includ- 
ing artificial flattening of the core density profile owing 
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Fig. 9. — A variety of disk-like KDCs at 15-20% gas. Top row: lR e velocity maps. Bottom row: Corresponding h-j maps. The 
color scaling and labels are as in previous figures. Each column is labeled at the top with the merger orbit and gas fraction, and at the 
bottom with the LOS. 



to the gravitational softening, and imprecise centering in 
the orbital classification, so we chose to omit the very 
innermost scales from the structural profiles. 

The remnants do not typically show large orientation 
twists within the 50th percentile in binding energy, or 
~ 1.5i? e (excluding the KDCs, since the shape may not 
be well-resolved on such small scales). The low f gas 
remnants typically have large kinematic misalignments, 
while the angular momentum of the 40% gas remnants 
is well aligned with the short axis of the potential. 

Profiles of the orbital structure, shape, and orientation 
of all 48 dissipative remnants (5-40% gas) are shown in 
Figures [15] to [20] of the Appendix. We have outlined the 
main trends with f gas in this section, but there is sub- 
stantial variation in the structural transformations from 
one merger orbit to the next. Some of the remnants (I 
and p) are already transformed into oblate, z— tube dom- 
inated systems at 5% gas, while others (e.g. i and k) are 
still box-dominated and substantially prolate at this gas 
fraction. Remnant o is oblate at 0% gas and becomes pro- 
late at 5% gas. A small gas component transforms rem- 
nant j from a nearly pure prolate spheroid dominated by 
a;— tube orbits, to a triaxial system with a large z— tube 
population. As the CMC alters the shape of the poten- 
tial, both directly and indirectly through its interactions 
with the stars, the phase space boundaries between the 
orbital classes move. This shifting of the orbital bound- 
aries can induce further changes in the potential, leading 
to a cascade effect. The causes of the more complex or- 
bital transformations in some of the individual remnants 
are an interesting topic for future study. 

3.3. Kinematically distinct cores and embedded disks 

The 15-40% gas remnants often contain a cold z— tube 
population that appears as a KDC or embedded disk in 
~ lR e velocity maps, similar to the sub -components ob- 
serve d in some of the SAURON galaxies (|Krainovic et al.l 
l2008h . 

This component typically appears as a KDC in the 15- 
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Fig. 10. — Orbital structure of a typical KDC remnant. 

Upper left: Velocity map of remnant i, 20% gas, viewed in an 
oblique projection. Upper right: LOSVD at location A. near the 
major axis within the KDC. Lower left: LOSVD at location B, near 
the major axis outside the KDC. Lower right: LOSVD at location 
C, near the minor axis. Color scaling and labels are as in previous 
figures. 



20% gas remnants, with a large kinematic twist relative 
to the stars around lR e . Figure [H] depicts five selected 
KDC systems, from a variety of different remnants and 
viewing angles. lR e velocity maps are presented in the 
top row, while the bottom row shows the correspond- 
ing maps of I13 to establish the disk-like character of the 
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KDCs. In some cases there are also large photometric 
twists on the KDC scale, as in column 4 of Figure [9l 
The KDCs are visible from nearly all viewing angles in 
the ~ 50% of the remnants that have them (see Fig- 
ures H3 - HU) ■ 

The outer parts of the KDC remnants show a variety of 
different types of kinematics, depending of the remnant 
and projection - some show no rotation (e.g. column 1 of 
Figure HJ, some minor-axis or oblique rotation (columns 
2 and 3), and some a significant amount of major-axis 
rotation either aligned with or counter to the KDC (but 
always less disk-like than the KDC; e.g. columns 2 and 
4). I13 and v are always tightly anticorrelated within 
the KDC, but may have a positive (columns 2 and 5) or 
weaker negative (column 4) correlation outside the core. 
Some of the remnants show multiple kinematically dis- 
tinct parts, e.g. remnant i at 20% gas displays minor-axis 
rotation with /13 and v correlated, and major-axis rota- 
tion counter to the KDC, with /13 and v anticorrelated, 
at the outer edge of the map. 

The cores of the KDC remnants are not nearly as dis- 
tinct in their underlyi ng orbital structure as they are 
in their kinematics, as Ivan den Bosch et all (|2008f ) ob- 
served in the dynamical modeling of NGC4365. This is 
because the mass of the streaming population produc- 
ing the KDC is generally small, and the KDCs consist of 
an oblate z— tube population on top of a system already 
oblate and dominated by z— tube orbits because of the 
strong CMC. Figure ITUl illustrates this point by showing 
the orbital breakdown of the LOSVDs at three different 
locations in the KDC remnant depicted in the second 

v (km / s) 



-max max 




-400 -200 200 400 
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Fig. 11. — Embedded disk formed by dissipation. Upper 
left: Velocity map of remnant i, 40% ga s, viewed at (9, <j>) = 
(72,77)°. Upper right: Same map with the new stars, formed fro 
m gas during the merger, removed. Lower left: Same map with 
only the 15% of the new stars that form ed last removed. Lower 
right: LOSVD at location A, broken down by orbital class and 
stellar ag e. Color scaling and labels are as in previous figures. 



column of Figure |H1 one within the KDC (point A) , one 
along the major axis outside the KDC, in a region show- 
ing no rotation (point B), and one near lR e along the 
minor axis, in the part dominated by minor-axis rota- 
tion (point C). The orbital population at all three loca- 
tions is (perhaps surprisingly) comprised of roughly the 
same mixture of z— tube, x— tube, and box orbits, with 
z— tube orbits always in the majority. The very differ- 
ent appearance of the three regions in the kinematic map 
owes to large differences in the degree of orbital stream- 
ing - at point A the z— tubes are highly streaming; the 
a:— tubes are highly streaming at point C; and at point B 
both classes of tubes are distributed symmetrically about 
v = 0. 

In the 30-40% gas remnants, the scale of the disk com- 
ponent is typically large enough to fill most or all of the 
SAURON-scale kinematic maps. The embedded disk in 
these remnants is typically quite cold, because the disk 
stars form in a slow trickle extending over ~2 Gyrs af- 
ter the merger ([Robertson et al.ll2006bl ; iHoffmanl l200l 
and are therefore subject to very little violent relaxation. 
Figure [TT] shows remnant i at 40% gas, which has one of 
the strongest embedded disks in our sample, at a view- 
ing angle about 30° from edge on. The contribution 
from old stars, present as collisionless particles in the 
initial disks, is shown in the upper right panel. When 
only old stars are included, the maximum velocity on 
the map drops from 135 to 66 km/s, and the rotation 
is no longer well-aligned with the major axis. The old 
stars do still display more major-axis rotation than most 
remnants at low f gas , since they also contract and spin 
up in response to the rapid condensation of the gas (e.g. 
Eggen et aLlll962t iBlumenthal et aTlll986t iGnedin et al.l 
2004t IVillalobos et al.H2009D . 

The lower right panel shows what is left when only the 
15% of the new stars that formed latest (~5% of the total 
stellar mass) are removed. The maximum velocity still 
plummets to 72 km/s, and what remains appears similar 
to a KDC system, since the most recently-formed stars 
dominate the outer portions of the embedded disk. The 
contributions of old and recently-formed stars to a typi- 
cal LOSVD near lR e in the embedded disk is shown in 
the lower right panel. The stars that formed late have 
a much smaller velocity dispersion, and when superim- 
posed on the hot old population they yield a velocity 
distributi on strongly skewed o pposite the mean velocity 
(see also Hoff man et al.ll2009aD . 

3.4. Outer structure 

We have shown that the outer parts of the remnants 
(outside ~ 1.5i? e ) are less relaxed than the inner parts, 
and arc largely unaffected by the gas content of the disks 
(see Figures and . While kinematic maps of gas- rich 
remnants are dominated by dissipational rotation on lR e 
scales, their kinematics farther out reflects the "dry" part 
of the merger remnant that is left essentially untouched 
by the gas, resulting in large kinematic twists between 
around 1 and 3i? e . 

This point is illustrated clearly in Figure rTJl which 
compares velocity maps of remnant m at and 40% 
gas on three different scales, 1, 3, and 7R e , in the same 
oblique projection. Profiles of the orbital populations are 
reproduced in the right-hand panel for comparison. On 
lR e scales the dissipationless and gas-rich remnants look 
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Fig. 12. — Kinematic transitions in the outer parts. The top row shows dissipationlcss remnant m, vie wed at (9,0) = (68,35)°, 
while the bottom row shows the same remnant at 40% gas, along th e same LOS. Left column: lR e velocity maps. Left center: 3R e velocity 
maps, with th e lR e maps superimposed inside the boxes. Right center: 7R e velocity maps, with the 1 and 3R e maps inside the boxes. 
Right: Radial profiles of the orbital content. All maps are plotted on the same velocity scale. Gaussian (rather than GH) fitting was used 
in this figure because of the nois y LOSVDs at large radii. 



entirely different, the former being a minor-axis rotator 
while the latter is dominated by rapid disk-like rotation. 
However when we zoom out to 3R e , the similarity in 
their outer structure becomes apparent. The gas has 
thoroughly transformed the inner part of the 40% gas 
remnant, while the outer parts are left as they were in 
the dissipationless case. The same observation can be 
made from the orbital profiles - the orbital populations 
outside ~ 1.5i? e are similar in the and 40% gas rem- 
nants, but are radically different within lR e , where the 
dissipationless remnant is comprised mostly of box or- 
bits and the 40% gas remnant is dominated by z— tubes. 
The similarity in the outer kinematic maps is even more 
apparent on 7R e scales. 

Sharp kinematic and orbital transitions such as those 
in Figure Q21 will be clearly observable in new sur- 
veys such as SMEAGOL, probin g the stellar o rbital 
structure out to ~ 3i? e scales (|Proctor et al.l 120091 : 
iFoster et all 120091 ) . or the PN.S survey, probing the 
kinematics out to ~ 7R e . scales us i ng planetary nebu- 
lae a s tracers (iDouelas et all 120021 : iRomanowskv et all 
200.4 ICoccato et alJl200fljT Predictions for the statistics 
of kinematic tw ists and abrupt chang es in the rotation 
parameter, Xr (|Emsellem et al.l I2007T ). i n such surveys 
will be further discussed and quantified in lHoffman et al.l 
(20093) • 

The intrinsic outer structure of the remnants is sum- 
marized in Figure I13[ which is the same as Figure [8] but 
for the structure between ~ 3 and 7R e - The orbital 
populations, shapes, and orbital streaming are far less 
sensitive to f gas than in the inner parts. The outer stel- 



lar population is dominated by tube orbits, as expected 
given the steep effective density profile at these radii, 
with a roughly even mix of x— tubes and z— tubes. The 
exchange between box and z— tube orbits with increasing 
fgas is still noticeable, though diminished. From visual 
inspection of some randomly selected orbits, box orbits 
in the halo tend to be more stochastic and spherical in 
shape (|Clozell 12008). and we will further investigate the 
nature of the outer box orbits in future work. 

Both classes of tubes have a higher rotation bias than 
in the inner parts. The a;— tube orbits are especially 
highly streaming, which is a signature of the stars' dy- 
namically cold origin. To highlight this point, we have 
also plotted the streaming fraction of the x— tube and 
z— tube orbits over the same range of relative binding 
energies for a series of re-mergers of the 20 and 40% gas 
remnants (see Hoff man et aTll2009al) . intended to repre- 
sent "dry" mergers between elliptical galaxies that have 
alread y exhausted their g a s in a previous major merger 
(e.g. Ivan Dokkuml 120051: iBell et al.1 120061 iNaab et all 
I2006bb iKormendv et all 120091 ). The streaming of the 
z— tubes is not that much different in the remnants of 
mergers between dynamically hot ellipticals than in the 
disk mergers, but the x— tube orbits show far less or- 
dered rotation in the dry merger remnants. In a merger 
between hot systems, the random motions of the stars 
provide a source of initial angular momentum about the 
long axis, with no preferred sense of rotation. The uni- 
form streaming of the outer x— tube orbits might be a 
tell-tale sign of a cold (late-type) progenitor in dynami- 
cal models of observed systems. 
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Fig. 13. — Outer structure (~ 3 — 7R e ), averaged over eight merger orbits. Same as Figure[8] except that the orbital structure 
is averaged over the 75-95th percentiles in binding energy, and the shapes and kinematic misalignment are evaluated within the 80th 
percentile. The intrinsic orientation twist is computed as max^j), with i and j running from the 25th to the 85th percentile in binding 
energy (~ 0.5 — 5R e ). On the tube rotation plot (upper right panel), we have added data points for re-mergers of the 20 and 40% gas 
remnants (labeled as "dry" mergers - see text). 



The shapes of the matter distribution between 3 and 
7R e are round (e ~ 0.8 on average), and nearly max- 
imally triaxial. The stars are substantially more flat- 
tened (e ~ 0.75) and oblate (T ~ 0.4) than the matter 
as a whole (including DM). Large kinematic misalign- 
ments are common in the outer parts, producing large 
kinematic twists between ~1 and 3i? e in the most gas- 
rich remnants, as shown previously in Figure 1121 The 
30-40% gas remnants also typically have large instrin- 
sic orientation twists between ~1 and 5i? e , owing to the 
rapid figure rotation of the inner component. 

4. SUMMARY AND CONCLUSIONS 

We have shown that a variety of observed kinematic 
structures can be accounted for just by varying the 
gas fraction in binary 1:1 mergers. The projected 
kinematics within lR e tends to fall into one of four 
categories: (i) rapid disk-like rotation about the major 
axis; (ii) a prominent, disk-like KDC with slow or 
minor-axis rotation farther out; (hi) uniform slow 
rotation, with v max < 40 km/s; and (iv) prominent 
minor axis rotation, with little or no rotation about 
the major axis. In Figure [14] we classify the eight 
remnants at each f gas into these four categories 



based on visual inspection of their velocity maps in 
projection along the y— axis (intrinsic variations far 
outweigh projection effects in the global appearance 
of the kinematic maps). The 0-10% gas remnants are 
generally slowly rotating or dominated by minor-axis 



rotation, in agreem e nt with the p reviou s 



results of 
iCox et all 



iNaab et all (|2006aD : Uesseit et ail (|2007t ) 
(2006^ The 15-20% gas remnants often look similar 
to the lower- f gas remnants but with disk-like KDCs 
at their centers, and bear an interesting resemblance 
to observed galaxies such as NGC4365 (Statlc ret all 
120041: Ivan den Bosch et all 120081: iHoffman et al.ll2009'cF 



NGC5813, and NGC4458 (lEnTsellem et all l2004t 
ICappellari et alll2007t fvan den Boschll2008[ ). 

The 30-40% gas remnants generically sho w rapid, disk- 
like r o tation with v/a up to ~ 1 (see also IRix fc W hite 
19901 IS corza fc Benderl 119951: iRobertson et all 12006b: 

None of them display the ex- 
' 2) of the most rapid rota- 
sample, although other stud- 



Krainovic et all 120081 ). 



treme rotation (v/a 
tors in the SAURO N 

ies (|Robertson et alj|2006bt IRobertson fc Bullock! 120081) 
have shown that 1:1 mergers between 60-80% gas disks 
can yield higher values of v/a. High v/a values are al- 
ways accompanied by high, anticorrelated /13 values in 
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our remnants, so the subset of SAURON rapid rotators 
with a shallow /13 — v j o relation may require a different 
formation mechanism. 

The variety of features that can be produced with 
1:1 mergers alone suggests that there are degeneracies 
in the formation scenarios leading to a given gross 
outcome. For instance cold embedded disks can arise 
either from gas-rich merging or fr om cosmological 
gas i nflow and secu l ar evolution (iNaab fc Burkertl 
20011: IRobertson et al.l l2004t iDekel fc Birnboiml 120061: 
Robertson et al.l 120061* iBournaud et al.l [2007b; 



Genzel et all 120081 : [Shapiro et al.l [20081: IDekel et al l 
20091: INaab et all 120091* IKeres fc Hernauistl 12001 
Ceverino et al. I 1200911 Lower progenitor mass ra- 
tios (~ 3-4:1 instead of 1:1 mergers) can produce 
remnants resembling many SAURON rapid rota- 
tors e ven in dissipa ti onless simulations (INaab et all 



199 a iCretton et al.l l200lt INaab fc Burkertl 



Bournaud et al. 2004, 2005b; Burkcrt & Naab 



2003 



2005 



Jesseit et all 120071: Hopkins et al.l 120091* Uesseit et all 
200a Uohansson et all 120091 ) . Sequential minor mergers 



can produce slowly rotating ellipticals as well as major 
mergers with low f gas , that are in better agreement 
with observed round, isotropic, and featureless sys- 
tems such as NGC4486 NGC4552, and NGC 5846 
dWeil fc Hernauistl [1991 [19951: IBournaud et al.l I2007at 
Cappellari et al.l 120071: IBurkert et all 120081: INaab et all 
2009a; iGonzalez-Garcia et al.l [2009). These scenarios 



must be distinguished w ith more detailed structural and 
dyna mi cal analysis (e. g. INaab et alj|2006bt Uesseit et all 
20071; IBurkert et all l2008t Ide Lorenzi et al.l 1 2008 : 



van den Bosch et al l 2008t IRobertson fc Bullockl I2OO8 
q^qt-,;™ of al I lonns* To^nn- ot al I lonnsl- iTrnffmc, ot d 



Shapi ro et all l2008t Ueong et al.l 12001 Hoffman et all 
2009a: iThomas et al.1120091: iRomanowskv et all 12009b; 



Romano-Diaz et al.l l2009f ). 



Our orbital analysis suggests a rather simple picture 
of how the structure of the remnants arises. Non- 
axisymmetric torques in the merger trigger the forma- 
tion of a triaxial structure with a large population of 
box orbits. When gas is present, the strong CMC fo rmed 
through dissipation (e.g. iMihos fc Hernqu ist 1994aj) con- 
verts the majority of the box population into z— tube or- 
bits, which have nearly canceling rotation since box dif- 
fusion has no preferred direction. Stars in inclined disks 
begin with a large polar angular momentum component 
(je) that is retained throughout the merger, so they do 
not cross orbital boundaries. This produces a streaming 
a:— tube population in the remnants, that is most pro- 
nounced at large radii. This picture of the x— tube pop- 
ulation is supported by the facts that (i) co-planar merg- 
ers produce remnants almost entirely devoid of a:— tube 
orbits; (ii) the phase space available to x— tube s in the 
remnant potentials i s generally underpopulated ([Barnes! 
Il992t iHoffrnanl 120071) : and (hi) in some observed ellipti- 
cals the velocity ellips oid is flattened along the polar axis 
van den Boscbll2008ri . 

In gas-rich merger remnants, there is also a pronounced 
streaming z— tube population on small scales owing to 
the gas that dissipates energ y and retains its angu lar mo- 
mentum during the merger ((Hopkins et aTll2009bl) . Out- 
side a fairly well-defined boundary (see e.g. Figures [7] 
and[T2j). the orbital structure is largely unaffected by the 
gas. This boundary corresponds to the radius where box 
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Fig. 14. — Classification of the remnant kinematics. The 

eight remnants of each gas fraction were categorized as (i) showing 
rapid major-axis rotation, (ii) slowly rotating with a disk-like KDC, 
(iii) slowly rotating throughout, or (iv) showing rapid minor-axis 
rotation, based on visual inspections of their velocity maps (see 
text). 

orbits cease to dominate in the dissipationless remnants, 
implying that these regions may be left relatively intact 
because they have never had a large population of radial 
orbits to convey information about the galactic center to 
their location. Time-dependent orbital classification and 
idealized evolving models are needed to verify this in- 
terpr etation (e.g. iHoffman et al|[2009cHVillalobos et all 
120091 ). 

The characteristic orbital structure arising from this 
picture, and in particular the sharp kinematic transitions 
predicted between ~ 1 and 3i? e , should be readily observ- 
able with surveys in progress such as SMEAGOL. Some 
pronounced kinematic transitions have been observed 
around these radii in surv ey pilot studies ([Proctor et all 
l2009llCoccato et alJl2009h . but more data will be needed 
to determine whether these transitions are of the same 
nature as those in the merger simulations. Note that 
some features, e.g. the uniform streaming of the a;— tube 
orbits, can only be captured with dynamical modeling - 
substantial minor-axis rotation alone could be produced 
by either a dominant population of a;— tubes with a small 
rotation bias, or a smaller, highly streaming population. 

We do not necessarily expect most real galaxies to dis- 
play the sharp features in their orbital structure found 
in this paper, since cosmological galaxy formation his- 
tories are far more complex than isolated 1:1 mergers 
between pure disks. However the results of our analy- 
sis illustrate how orbital analysis can isolate subsets of 
the stellar population with similar histories and place in- 
tuitive constraints on galaxy formation mechanisms. In 
future work we hope to extend this type of analysis to a 
broader range of merger parameters and more complex 
cosmological formation scenarios. 
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APPENDIX 

In Figures [T5l - |2T)1 we present radial profiles of the orbital populations, intrinsic shapes, and orientations of all eight 
merger remnants with each initial gas fraction (5, 10, 15, 20, 30, and 40%). The format of the figures is identical to 
that of Figure U] in the body of the paper (instrinsic structure of the dissipationless remnants). The figures may be 
viewed one-at-a-time to get a feel for the variation with merger orbit at fixed f gas , or compared row-for-row to follow 
the evolution of the structure with f gas for a fixed merger orbit . 
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Fig. 15. — Intrinsic structure of the 5% gas remnants. 
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Fig. 16. — Intrinsic structure of the 10% gas remnants. 
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Fig. 17. — Intrinsic structure of the 15% gas remnants. 
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Fig. 18. — Intrinsic structure of the 20% gas 



remnants. 
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Fig. 19. — Intrinsic structure of the 30% gas remnants. 
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Fig. 20. — Intrinsic structure of the 40% gas remnants. 



